Code
library(zoib)
library(tidyverse)
library(GGally)
library(ggpubr)
library(kableExtra)
library(mice)
library(parallel)
theme_set(theme_pubr(legend = "bottom"))library(zoib)
library(tidyverse)
library(GGally)
library(ggpubr)
library(kableExtra)
library(mice)
library(parallel)
theme_set(theme_pubr(legend = "bottom"))nnns_imputed<- readRDS("../Haojia_work/nnns_imputed.rds")
dat <- lapply(1:20, function(i) complete(nnns_imputed, i))\[ f(y_i|\eta_i) = \begin{cases} p_i, & y_i =0 \\ (1-p_i)q_i, & y_i =1 \\ (1 − p_i)(1 − q_i)Beta(\alpha_{i_1}, \alpha_{i_2}) &y_i \in (0, 1) \end{cases} \]
\[ \begin{aligned} logit\left(\mu^{(0,1)} = \frac{\alpha_1}{\alpha_1+ \alpha_2}\right) = \mathbf x_1 \boldsymbol \beta_1 + I_1(\mathbf z_{1,i}\gamma_{1,i}) \\ log\left(v_i = \alpha_1+ \alpha_2\right) = \mathbf x_2 \boldsymbol \beta_2 + I_2(\mathbf z_{2,i}\gamma_{2,i}) \\ logit\left(p_i\right) = \mathbf x_3 \boldsymbol \beta_3 + I_1(\mathbf z_{3,i}\gamma_{3,i}) \\ logit\left(q_i\right) = \mathbf x_4 \boldsymbol \beta_4 + I_1(\mathbf z_{4,i}\gamma_{4 ,i}) \end{aligned} \]
set.seed(11282023)
tictoc::tic()
# Define a function for model fitting
fit_model <- function(d) {
zoib(
Percent_of_feeds_taken_by_mouth_at_discharge ~
Pre_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female #x1 design matrix
| 1 | #x2 design matrix
Pre_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female | #X3 design matrix
Pre_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female, #x4 design matrix
data = d,
n.response = 1,
zero.inflation = TRUE,
one.inflation = TRUE,
link.mu = "logit",
link.x0 = "logit",
link.x1 = "logit",
random = 0,
n.chain = 4,
n.iter = 3000,
n.thin = 2,
n.burn = 200,
seeds = c(11, 29, 20, 23)
)
}
model_results <- lapply(dat, fit_model)
# Save results
saveRDS(model_results, "pre_op_models.rds")
tictoc::toc()set.seed(11282023)
tictoc::tic()
# Define a function for model fitting
fit_model <- function(d) {
zoib(
Percent_of_feeds_taken_by_mouth_at_discharge ~
Post_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female #x1 design matrix
| 1 | #x2 design matrix
Post_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female | #X3 design matrix
Post_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female, #x4 design matrix
data = d,
n.response = 1,
zero.inflation = TRUE,
one.inflation = TRUE,
link.mu = "logit",
link.x0 = "logit",
link.x1 = "logit",
random = 0,
n.chain = 4,
n.iter = 3000,
n.thin = 2,
n.burn = 200,
seeds = c(11, 29, 20, 23)
)
}
# Parallelize model fitting
model_results <- lapply(dat, fit_model)
# Save results
saveRDS(model_results, "post_op_models.rds")
tictoc::toc()b: vector of estimates from Eqn 1; that is, g(mu) = xb*b +z*gamma
d: vector of estimates from Eqn 2; that is, log(eta) = xd*d+z*gamma
b0: vector of estimates from Eqn 3; that is, g(p0) = x0*b0 +z*gamma
b1: vector of estimates from Eqn ;4 that is, g(p1) = x1*b1+z*gamma
post_op_models =readRDS(file="post_op_models.rds")
pre_op_models =readRDS(file="pre_op_models.rds")
post_op_coeff = list()
pre_op_coeff = list()
for(i in 1:length(post_op_models)){
pre_op_coeff[[i]] = pre_op_models[[i]]$coeff
post_op_coeff[[i]] = post_op_models[[i]]$coeff
}
pooled_pre_op = runjags::combine.mcmc(mcmc.objects = pre_op_coeff, collapse.chains = FALSE)
pooled_post_op = runjags::combine.mcmc(mcmc.objects = post_op_coeff, collapse.chains = FALSE)
saveRDS(pooled_pre_op, "pooled_pre_op.rds")
saveRDS(pooled_post_op, "pooled_post_op.rds")pooled_pre_op = readRDS("pooled_pre_op.rds")
pooled_post_op = readRDS("pooled_post_op.rds")pooled_pre_op |> traceplot()pooled_post_op |> traceplot()pooled_pre_op |> autocorr.plot()pooled_post_op |> autocorr.plot()summary_pooled_pre_op = pooled_pre_op |> summary()
summary_pooled_post_op = pooled_post_op |> summary()
rnames = c("Baseline","Attention Score","Length of Intubation (d)",
"Single Ventricle w/ Arch Obstruction","Two Ventricles w/ Arch Obstruction",
"Age", "Female")
pre_mean = summary_pooled_pre_op$statistics[,"Mean"]
pre_lb = summary_pooled_pre_op$quantiles[,"2.5%"]
pre_ub = summary_pooled_pre_op$quantiles[,"97.5%"]
post_mean = summary_pooled_post_op$statistics[,"Mean"]
post_lb = summary_pooled_post_op$quantiles[,"2.5%"]
post_ub = summary_pooled_post_op$quantiles[,"97.5%"]
pre_b_df= cbind("Mean"= pre_mean[1:7] |> exp() |> round(2),
"2.5%"= pre_lb[1:7] |> exp() |> round(2),
"97.5%" = pre_ub[1:7] |> exp() |> round(2))
pre_b0_df= cbind("Mean"= pre_mean[8:14] |> exp() |> round(2),
"2.5%"= pre_lb[8:14] |> exp() |> round(2),
"97.5%" = pre_ub[8:14] |> exp() |> round(2))
pre_b1_df= cbind("Mean"= pre_mean[15:21] |> exp() |> round(2),
"2.5%"= pre_lb[15:21] |> exp() |> round(2),
"97.5%" = pre_ub[15:21] |> exp() |> round(2))
pre_df = cbind(pre_b_df,pre_b0_df,pre_b1_df)
rownames(pre_df) = rnames
post_b_df= cbind("Mean"= post_mean[1:7] |> exp() |> round(2),
"2.5%"= post_lb[1:7] |> exp() |> round(2),
"97.5%" = post_ub[1:7] |> exp() |> round(2))
post_b0_df= cbind("Mean"= post_mean[8:14] |> exp() |> round(2),
"2.5%"= post_lb[8:14] |> exp() |> round(2),
"97.5%" = post_ub[8:14] |> exp() |> round(2))
post_b1_df= cbind("Mean "= post_mean[15:21] |> exp() |> round(2),
"2.5%"= post_lb[15:21] |> exp() |> round(2),
"97.5%" = post_ub[15:21] |> exp() |> round(2))
post_df = cbind(post_b_df,post_b0_df,post_b1_df)
rownames(post_df) = rnames
pre_kable = pre_df |>
kable(caption = "Odds Ratio of Percent Oral Feed Model Results (Pre-Operation)") |>
add_header_above(header = c("Predictor" = 1,
"Odds of Oral Feed \n when oral feed between 0 and 1" = 3,
"Odds of 0% Oral Feed" = 3,
"Odds of 100% Oral Feed" =3)) |>
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) |>
add_footnote(paste("Posterior variance is estimated to be ", pre_mean[22] |> exp() |> round(2),
" (",pre_lb[22] |> exp() |> round(2),",", pre_ub[22]|> exp()|> round(2),")"))
post_kable = post_df |>
kable(caption = "Odds Ratio of Percent Oral Feed Model Results (Post-Operation)") |>
add_header_above(header = c("Predictor" = 1,
"Odds of Oral Feed \n when oral feed between 0 and 1" = 3,
"Odds of 0% Oral Feed" = 3,
"Odds of 100% Oral Feed" =3)) |>
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))|>
add_footnote(paste("Posterior variance is estimated to be ", post_mean[22] |> exp() |> round(2),
" (",post_lb[22] |> exp() |> round(2),",", post_ub[22]|> exp()|> round(2),")"))pre_kable| Mean | 2.5% | 97.5 | Mean | 2.5% | 97.5 | Mean | 2.5% | 97.5 | |
|---|---|---|---|---|---|---|---|---|---|
| Baseline | 0.96 | 0.21 | 4.47 | 0.20 | 0.02 | 2.32 | 0.02 | 0.00 | 1.11 |
| Attention Score | 1.02 | 0.71 | 1.47 | 0.74 | 0.37 | 1.35 | 1.74 | 0.63 | 5.34 |
| Length of Intubation (d) | 0.86 | 0.76 | 0.98 | 1.35 | 1.13 | 1.67 | 0.69 | 0.45 | 0.99 |
| Single Ventricle w/ Arch Obstruction | 1.57 | 0.71 | 3.36 | 6.19 | 1.90 | 22.14 | 0.76 | 0.03 | 10.62 |
| Two Ventricles w/ Arch Obstruction | 0.81 | 0.42 | 1.59 | 3.43 | 1.11 | 11.58 | 0.52 | 0.08 | 3.05 |
| Age | 0.98 | 0.92 | 1.04 | 0.97 | 0.87 | 1.07 | 1.19 | 1.02 | 1.43 |
| Female | 0.65 | 0.37 | 1.13 | 0.56 | 0.21 | 1.44 | 1.94 | 0.40 | 9.96 |
post_kable| Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | |
|---|---|---|---|---|---|---|---|---|---|
| Baseline | 0.39 | 0.04 | 3.29 | 0.03 | 0.00 | 0.67 | 0.02 | 0.00 | 37.09 |
| Attention Score | 1.20 | 0.82 | 1.75 | 1.19 | 0.71 | 2.01 | 1.33 | 0.34 | 4.09 |
| Length of Intubation (d) | 0.88 | 0.77 | 1.00 | 1.34 | 1.12 | 1.65 | 0.72 | 0.47 | 1.03 |
| Single Ventricle w/ Arch Obstruction | 1.72 | 0.77 | 3.71 | 6.20 | 1.94 | 21.52 | 0.84 | 0.02 | 13.07 |
| Two Ventricles w/ Arch Obstruction | 0.85 | 0.44 | 1.64 | 3.16 | 1.06 | 9.91 | 0.63 | 0.09 | 3.92 |
| Age | 0.98 | 0.92 | 1.05 | 0.97 | 0.87 | 1.07 | 1.21 | 1.02 | 1.45 |
| Female | 0.64 | 0.36 | 1.11 | 0.54 | 0.20 | 1.40 | 2.47 | 0.46 | 14.76 |